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Using numerical and analytical methods implemented for different models we conduct a systematic 
study of thermodynamic properties of pairing correlation in mesoscopic nuclear systems. Various 
quantities are calculated and analyzed using the exact solution of pairing. An in-depth comparison 
of canonical, grand canonical, and microcanonical ensemble is conducted. The nature of the pairing 
phase transition in a small system is of a particular interest. We discuss the onset of discontinuity 
in the thermodynamic variables, fluctuations, and evolution of zeros of the canonical and grand 
canonical partition functions in the complex plane. The behavior of the Invariant Correlational 
Entropy is also studied in the transitional region of interest. The change in the character of the phase 
transition due to the presence of magnetic field is discussed along with studies of superconducting 
thermodynamics. 
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I. INTRODUCTION 

Pairing correlations and related superconducting or su- 
perfluid properties are robust features of quantum many- 
body systems. In physics anywhere from quarks to stars 
it is hard to find systems that under certain conditions 
do not exhibit pairing correlations. The Cooper phe- 
nomenon [3], namely the instability against formation of 
particle-pairs in a macroscopic Fermi-system under an 
arbitrarily weak attractive force, is a primary reason for 
thriving of pairing. 

Pairing in mesoscopic systems, such as atomic nuclei 
metal clusters 0,0; 01, ultra small grains quan- 
tum dots 0, interacting spins [!,[§], has attracted a lot 
of attention rece ntly . Indeed, questions of phase tran- 
sitions ID, floL [Til . fl2l. fT3|. interplay with other collective 
modes continuum effects and thermodynam- 

ical properties of small systems are important for the 
| present day science and technology. 

In this work we conduct a systematic study of thermo- 
dynamics of pairing correlations in small systems. We 
use two-types of model Hamiltonians of lower and higher 
symmetry where the pairing problem is solved exactly 
and all quantum states are identified. We use a quasi- 
spin algebra with the effective numerical implementation 
to obtain a full solution for systems ranging in size from 
a few particles to as large as over a hundred of parti- 
cles. The traditional BCS solution is also considered 
for comparison. Using these results we compare differ- 
ent thermodynamic ensembles: microcanonical, canoni- 
cal and grand canonical. The differences indicate a meso- 
scopic nature of the system [l7], Ell EH] and diminish in 
the macroscopic limit. Some discrepancies observed in 
thermodynamics are related to non-thermal nature of the 
pure pairing interaction [20] and raise questions of equi- 
libration and thermalization. Through thermodynamic 
ensembles and using invariant correlational entropy we 
study and analyze the pairing phase transition as a func- 
tion of temperature or excitation energy, magnetic field, 
size of the system, and pairing strength. 



We further explore the evolution of zeros in the 
complex temperature plane for the canonical ensemble 
[2ll I22I I23L [24 |. where recent findings established clear 
correlations of pair breaking with peaks in entropy and 
branches of complex temperature roots approaching real 
axis @, H, We extend this discussion with consid- 
eration of the phase transition based on the Yang-Lee 
theory [IE H3, El| • The study of the system in the mag- 
netic field, evolution of zeros in the partition function 
as a function of the field strength, spin fluctuations and 
the change of the phase transition type are particularly 
interesting. 

The presentation below is structured as follows: we 
first introduce the pairing Hamiltonian, identify prop- 
erties of the pairing problem and define models for our 
study in Sec. Ull In Sec. [TTT] we consider a BCS ap- 
proximation which shows the generic features of a paired 
system. The bulk of the work is presented in Sec. HVI and 
its subsections, where different methods are introduced, 
discussed and compared. In Sec. [Vl]we concentrate on 
the effects that external magnetic field or rotation have 
on the properties of paired systems; this includes the clas- 
sification of the phase transitions using the distribution 
of zeros in the partition functions. 



II. PAIRING HAMILTONIAN 

We approach the pairing problem by defining a pair of 
two single-particle states denoted here as 1 and 1. This 
pair-wise identification can be based on an arbitrary sym- 
metry; however, the fundamental symmetry with respect 
to time reversal is the most common. For this work we 
assume a pair as two particles in time-conjugated single- 
particle states that due to this symmetry have identi- 
cal energies. Using the language of the second quanti- 
zation the pair creation and annihilation operators are 
p\ = a\a- and p\ = ajai, respectively. Here the a\ and 
ai are single-particle creation and annihilation operators 
with the usual fermion commutation rules. The pair is 
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labeled by the same single-particle index 1, and is invari- 
ant under the time conjugation, p\ = pj, since a~ = a. 

The algebra of the pair operators on a pair-state 1 (a 
pair of orbitals 1 and 1) is identical to that of an SU(2) 
spin algebra called quasi-spin, in general the commuta- 
tion relations are 

\pt,p 2 ]=26 12 pi, (1) 

where 

Pi = ("i - ^ , (2) 

the operator related to the particle number n\ — (a\a\ + 
a~aj)/2 operator for the pair-state 1. 

A pair-state (1,1) occupied by a pair or completely 
empty correspond to quasi-spin 1/2 with projections 
p\ = 1/2 and pf = —1/2, respectively. Alternatively, 
these states are referred to as states with seniority s\ = 
identifying the number of unpaired nucleons in the pair- 
state 1. The states with one unpaired particle correspond 
to si = 1 and to zero quasi-spin. 

The most general form of the two-body Hamiltonian 
that describes motion of pairs at fixed particle number is 

H = 2j2cmi- E G ^P\P2, (3) 

1>0 1,2>0 

where the summation runs over pair-orbitals, denoted as 
1 > 0, ei are single-particle energies, and Gn = G21 
determines the strength of pair scattering. Using the 
quasi-spin the same Hamiltonian can be written as 

H = E £ i+E 2 ( £ i - p z i- E G ^ (ft • ft - *m ■ 

1>0 1>0 ^ ' 12>0 

(4) 

The problem is analogous to the Heisenberg model of 
0/2 — s interacting spins \p\ = 1/2, with the Zeeman 
splitting created by the single-particle energies. The 0/2 
stands here for the total number of double-degenerate 
levels and s — J2i s i represents the total seniority. Be- 
cause of the magnetic-field like splitting the total quasi- 
spin vector p — Xa>oft ls not conserved, while the re- 
maining cylindrical symmetry allows for the conservation 
of the z-projection p z = N/2 — 0/4 equivalent to the to- 
tal particle number TV = 2^ 1>0 m. 

The eigenstates of the Hamiltonian J3]) are identified 
by the set of 0/2 seniorities s = {si} denoting the avail- 
able and blocked pair-states. In the language of the spin 
model ((4]) seniorities represent the number of spin 1/2 
particles in the system, thus totally removing all blocked 
states from interaction. The Hamiltonian within a cer- 
tain seniority partition s is given as 

i7 s = ^ Sl£l + 2E«i^i-^ i ) (5) 

1>0 1>0 ^ ' 

s 

~E Gl2 (ft 'Pz-PiPi) 1 

1/2 



where the upper summation limit s implies that all 
blocked states with si = 1 are excluded. 

Since each unpaired particle doubles the degeneracy 
of the many-body state the total degeneracy of a given 
eigenstate is g s — 2 s . With other symmetries, beyond 
the time reversal, the degeneracy of states can be higher. 
Additional degeneracies such as the one due to the rota- 
tional symmetry can further reduce the problem to larger 
values of the quasi-spin. In the spherical shell model 
within a given j-shell there are total of uij = j + 1/2 
time-conjugate pair states, and the total quasi-spin is 
preserved by the pairing interaction. For such j-shell 
a quasi-spin vector Pj = h^2 m Pjm CeLH be introduced 
which together with the number operator for this level 
and its own hermitian conjugate again forms an SU(2) 
group. The independence of matrix elements and quasi- 
spin operators on magnetic sub-states allows to rewrite 
the Hamiltonian ([3]) as 

// E f < v < E 1 - 7 ''' - (6) 

3 n' 

where for the reasons of the two-particle state normaliza- 
tion a pair operator and interaction matrix elements are 
redefined as follows 

V J m 

Vjj< = y/UJjtVj'Gjj' . 

The exact diagonalization of the pairing Hamiltonian ([3]) 
or (J6j) , depending on the symmetries of the model, is 
performed using the quasi-spin algebra. The ability to 
obtain all many-body states with a relatively simple ex- 
act treatment of pairing is an important component in 
this study. The more detailed discussion of the seni ority 
based diagonalization can be found in Ref. [1^, [3(J l3lf . 
We refer to the exact treatment of pairing as EP. The 
applications of algebraic methods extend far beyond our 
models; treatments of proton-neutron pairing as well as 
more exotic forms of p airing-type Hamiltonians are dis- 
cussed in Ref. [13, El H \M \M, S EH, EH 53] ■ Other 
methods of exact solution, analogies with boson-fermion 
models and electrostatic analogies should be mentioned 

S3, S3, S^i . 

Below in Sec. (|IVp we introduce thermodynamic en- 
sembles and discuss thermodynamic variables used to 
study the many-body system that undergoes pairing 
phase transition. For each of the cases we construct 
the partition function exactly based on the full numer- 
ical solution to the pairing problem. As our examples 
we consider two basic types of systems. The picket-fence 
(or ladder system) which has 0/2 equally spaced double- 
degenerate levels, where the total fermion capacity is O. 
The level spacing is chosen as the unit of energy. The 
picket-fence model is a minimal symmetry system with 
the time reversal only; therefore the degeneracy of each 
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eigenstate a is g as = 2 s . A second model with only two 
levels, but of large degeneracy, represents an opposite 
"high symmetry" case. Spacing between the two levels 
is again taken as the unit of energy. Due to additional 
symmetry, the degeneracy of many-body states is higher. 
The action of the pairing Hamiltonian is limited to either 
diagonal or level to level pair transfer. For the two-level 
system with an appropriate selection of the basis states 
the Hamiltonian matrix is tri-diagonal. This facilitates 
substantially the numerical treatment, making determi- 
nation of all many-body states in systems with a hundred 
or more particles possible. The two types of model spaces 
with total occupancy O, the particle number N, and the 
constant pairing strength G constitute the set of input 
parameters in this study. Introduction of the magnetic 
field in Sec lVIl does not require a separate diagonaliza- 
tion, however requires determination of the total spin 
projections onto an axis parallel to the direction of the 
field. We note that the total number of many-body states 
is o - 0! 



III. BCS 

The BCS approximation is the common approach to 
tackle the pairing problem. While this method is asymp- 
totically exact in thermodynamic limit it still produces 
remarkably good results for smaller systems. The BCS 
method assumes the presence of a condensate and ap- 
proximates the dynamics of interacting particles ([3]) with 
a motion of independent quasi-particles. Although most 
of the issues that we intend to address in this work can 
not be fully explored within the BCS picture due to its 
limitations, the method is a good benchmark for many of 
the questions and an excellent guidance to the dynamical 
regions of interest. Below we review the approach while 
stressing some of the key elements relevant to this work. 

Within the BCS theory the general pairing Hamilto- 
nian in Eq. ((3J) is brought to an approximate single par- 
ticle form using Bogoliubov's transformation. The pa- 
rameters of the transformation, the set of gaps Aiand 
chemical potential /Lt, are determined via gap equations 

A 2 

(8) 



The total energy of the paired system is 



2>0 



and the chemical potential is given by the particle num- 
ber 



N = 2^^7ii where ri\ = — (1 



1>0 



^ 

2 e\ 



(9) 



For simplicity of notations we introduce single particle 
energies shifted by the chemical potential and the diago- 
nal interaction strength e\ = e% — fi — Gu/2. The result of 
the Bogoliubov transformation is the spectrum of states 
given by the independent quasi-particle excitations with 
energies 



e x = Je\ + A? . 



(10) 



E : 



1>0 



Gn 
~~2~ 



ni 



1,2>0 



AiA 2 

> 

' 4eie 2 



As earlier, the summations here go over the pair-states. 

In this work for all our models we use a constant pair- 
ing strength Gu = G which due to Eq. |(8]) leads to a 
constant pairing gap for all single particle pairs, Ai = A. 
A single parameter for the interaction strength, in our 
view, allows for the most transparent study of the impor- 
tant features, the results are generic, and the methods of 
BCS and EP are applicable to general situations. For 
constant pairing the BCS gap equation and the energy 
are a textbook examples: 



2 ^ ei 

i>o 1 



E 



1>0 



G 



ni 



G ' 



(11) 



To accommodate the cases with higher symmetry follow- 
ing Eq. ([7]) it is convenient to introduce V — uoG where 
w is the pair degeneracy which is level independent in 
both picket-fence (u)j = 1) and two-level (ujj 1 = u)j 2 = uj) 
models. 

The particle number non-conservation intrinsic to the 
Bogoliubov transformation is one of the problems asso- 
ciated with the BCS applications to mesoscopic systems. 
Furthermore, in a system with discrete levels Eq. ([8]) 
may not have a solution, with the exception of a trivial 
case Ai = 0. Formally, this transitional point cor- 
responds to the critical interaction strength where the 
largest eigenvalue of the matrix built from the elements 
(Gi2£i+G2i£2)/(4ei£2) is equal to unity. The interpreta- 
tion of this is that at a low pairing strength the pairing is 
too weak to overcome gaps in the single particle spectrum 
which leads to a normal state. This situation is again spe- 
cific to small systems where it appears in contrast to the 
Cooper instability [l[. The total absence of the pairing 
correlations below the critical pairing strength is a sec- 
ond major drawback of the BCS approach in mesoscopic 
systems. Exact solutions indicate a gradual dissipation 
of pairing correlations extending almost to zero strength 
[3, H3, [4§, The critical pairing strength as deter- 

mined by the BCS is still an important parameter iden- 
tifying the location of the mesoscopic phase transition. 

An analytic solution to the BCS equations can be 
obtained for the system of two levels defined above. 
For a half-occupied system the chemical potential due 
to the particle-hole symmetry is an exact average of 
the monopole-renormalized single-particle energies /i = 
(ei + e 2 ~ G)/2. Thus, 



a- = r- - I — 



Ae = e\ - e 2 . 



(12) 



The introduction of the renormalized strength V makes 
this equation independent of O. 

In Fig. QJa) the BCS gap is plotted as a function of 
energy for a two- level model following Eq. (fl2]l . The 
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curve has a square-root discontinuity at the critical pair- 
ing strength ^=0.50 in the units of level spacing. The 
concept of the gap does not appear in the exact solution, 
however this quantity can be deduced from the energy 
associated with paring correlations. The second curve 
in the same figure shows the gap computed through Eq. 
(fTTjl where energy and occupation numbers are obtained 
from the exact solution. The difference between these two 
curves depicts the shortcoming of the BCS when applied 
to a small system; for related discussions and comparison 
of BCS with exact techniques see Refs. [13, E3, E2] ■ In 
Fig. \V[b) an alternative view on the EP-BCS comparison 
is given. Here we show the energy difference per particle 
between BCS and EP as a function of the pairing strength 
for N=20 and 100 particles. As the particle number grows 
the BCS and EP become equivalent. The peak in the 
BCS-EP discrepancy appears in the pairing phase tran- 
sition region, around V cr « 0.6 which is close to an ana- 
lytically obtained BCS value of 0.5. The discrepancy in 
V cr is known to arise from the pair-vibrations and other 
renormalizations of the BCS ground state [H, E3|- 

For our second (picket-fence) model, the critical pair- 
ing strength can be determined in the case of a half- 
occupied system with an even number of levels through 
the sum of a harmonic series 
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(13) 



which in the limit of a large number of levels converges 
to zero logarithmically G cr ~ Ae/lnfi. We remind here 
that formally for this model G = V since u>j — 1 and £7/2 
equals to the number of levels. This logarithmic depen- 
dence in the macroscopic limit is related to the exponen- 
tial dependence of the gap on the pairing strength and 
density of states near the Fermi surface, which represents 
the Cooper instability. 

We conclude this section with a note on the BCS ap- 
proach at finite temperature T = 1//3. By modeling the 
thermodynamics of quasiparticles with non-interacting 
Fermi gas we obtain a modified version of the Eq. (fTTjl 



1 = ?E 



tanh 



(«*) 



1>0 



ei 



(14) 



where quasiparticle energies are of the form l|10p . A re- 
lated discussion of thermodynamic treatment within the 
grand canonical partition function and following from it 
thermal BCS is presented in Sec. IIVBI 

For the two-level, half occupied model the temperature 
dependence of the critical pairing strength is given by 



V cr = — COth( — ). 



(15) 



Figure 1: (Color online) In the upper panel the BCS pairing 
gap is shown as a function of the pairing strength for the two- 
level, half-occupied system with 20 particles. In the lower 
panel the energy difference per particle between BCS and the 
exact result is shown as a function of the pairing strength for 
the same JV=20 system and is compared with the results for a 
larger half-occupied two-level model containing 100 particles. 



IV. STATISTICAL TREATMENT 

Statistical properties of a many-body systems are ad- 
dressed using normalized density operators [H3|, usually 
referred to as statistical operators w [56] and defined as 



w(E, N) — — S(E - H)S(N - N) 



for the microcanonical, 



1 



(13, N) = - exp(-PH)5(N-N) 



for the canonical, and 



>(/?,/i) = i exp(-l3{H-fiN) 



(16) 



(17) 



(18) 



for the grand canonical ensemble. In the above defini- 
tions the parameter f3 = 1/T refers to an inverse temper- 
ature and /i corresponds to the chemical potential. Here 
we use units where the Boltzmann constant is equal to 
unity, allowing units of energy to be used for tempera- 
ture. The normalization constants Z, Z, and Z are the 
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partition functions for the corresponding ensembles; so 
that the statistical operators are normalized by the trace 
Tr(w) — 1. The statistical averages are calculated as 

(6) = Tr(6w). (19) 

The entropy for the above ensembles is defined as 

S = -(ln(tu)) = -Tr(w Into) . (20) 

The above definition is strictly speaking applicable only 
for thermally equilibrated system which makes ther- 
modynamical Boltzmann-Gibbs entropy discussed below 
equivalent to the von-Neumann entropy of a quantum 
ensemble in Eq. lf20|) . A new light on complexity of 
quantum states in non-thermalized or non-equilibrated 
systems can be obtained with the invariant correlational 
entropy [HtJ (ICE) that also appears to be a good tool 
to study the phase transitions in mesoscopic systems 
[HsL [H§|. The correlational entropy is defined through 
the behavior of the microcanonical density matrix lfl6|) 
for each individual quantum state in response to a noise 
in an external parameter. For the purposes of this work 
we consider pairing strength V to be this external param- 
eter. The variations in V within the interval [V, V + SV] 
result in an averaged density operator 

^ pV+SV 
™«=gyj v ™ a(n 

where the weight operator w a is a density operator for 
an individual quantum state a followed with evolution of 
V, for a fixed parameter V this is a projection operator. 
The averaged statistical weight matrix is used to obtain 
the ICE via Eq. J20]). 

The quality or applicability of a given thermodynamic 
approach to a small system is often under question. 
While in some studies various ensembles are used inter- 
changeably, there are significant dangers on this path. 
Our investigations below not only show up the pairing 
phase transition and its evolution as a function of the 
particle number but also draw attention to some subtle 
differences in thermodynamic treatments. 



A. Canonical ensemble 

Given an exact solution to the pairing problem via di- 
agonalization in the seniority scheme, Sec. HU the formal 
definition 1(17)1 can be written explicitly for the eigen- 
states labeled by a and s: 

w as = — cxp(-fjE as ) , where (21) 



Z{(3, N)=J2 9as eM~PE a , s ) (22) 

as 



is the canonical partition function. The ensemble average 
(fl9)l for any quantity is given as 

(°) = ^29asW as (as\0\as), (23) 

as 

where (as|0|as) is the quantum-mechanical expectation 
value for the corresponding operator in the eigenstate a 
with the seniority set s. The entropy is given via the 
usual expression 

S = ~y] g as w as ln(w aa ). (24) 

CVS 

The reader may be familiar with the following set of tra- 
ditional thermodynamic relations [sq| 

(£} = -Aln(Z), (25) 

the entropy S can be found directly from the statistical 
definition {20} 

S = lnZ + 13(E) = ~. (26) 

The Helmholtz free energy is defined as 

F = -Tln(Z) = (E) - TS. (27) 

The Eq. {25)1 involves a derivative, however in our calcu- 
lations we avoid numerical differentiations always going 
back to the definition ((23)l . For example specific heat 
is computed using its relation to the energy fluctuations 

({E-(E)f), 

The results of our study based on the canonical ensem- 
ble are shown in Fig. [2][H1 In Fig. [2] (a - d) free energy, 
entropy, energy, and energy fluctuation of the ladder sys- 
tem with 12 levels and 12 particles are shown as a func- 
tion of temperature, similar study may be found in @,[6l| 
and references therein. The critical pairing strength for 
this model from BCS, Eq. (fT3)l . at zero temperature is 
V cr = 0.27. The curves correspond to different pairing 
strengths showing various conditions: weak pairing with 
about half the critical pairing strength V = 0.13; pair- 
ing strength above the critical value V = 0.6; and strong 
pairing V = 1. All of the plots show essentially similar 
trends: there is a sharp change in each of the quantities 
as a function of temperature in a certain region. This 
region is associated with the phase transition from the 
paired to the normal state. Most transparently it can 
be seen in [2je) where it is associated with the peak in 
heat capacity. The critical temperature T cr depends on 
the pairing strength. It can be observed that the tran- 
sitional region for strong pairing (V > V cr for T = 0) is 
roughly consistent with the BCS, which gives T cr = 2.7 
and 1.3 for V = 1 and 0.6, respectively. Naturally, the 
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stronger pairing interactions support the superconduct- 
ing state at higher temperature or excitation energy. For 
weak pairing the transitional behavior is present at zero 
temperature. This is consistent with the earlier finding 
that pairing correlations appear in the ground state even 
for small V. The decline of weak pairing (V < 0.13) 
phase is still associated with the peak in heat capacity 
which becomes smaller as the pairing strength is weak- 
ened, while staying essentially at the same T cr ~ 1.3. 
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Figure 2: (Color online) (a) Free energy, (b) Entropy, 
(c) Energy, (d) Energy fluctuations, (e) Specific heat, and (f) 
Order parameter of a ladder system with 12 levels and 12 
particles as a function of temperature. 

The phase diagram can be further explored by consid- 
ering an order parameter which we define here as a frac- 
tion of paired particles ip — (N — (s))/N, the (s) is the 
ensemble-averaged value of the total seniority. The de- 
pendence of the order parameter on temperature, shown 
in Fig. [2](f ) , shows that the fraction of superconducting 
pairs drops sharply in the transitional region which is 
also identified by the critical behavior of other thermo- 
dynamic quantities. 

The contour plot of the order parameter as a function 
of the pairing strength and temperature is shown in Fig. 
[3j The shaded area in the upper left corresponds to the 
high percentage of particles in the condensate, which oc- 
curs at low temperature and high pairing strength; while 



in the opposite limit the superconducting state disap- 
pears. The solid line indicates the phase boundary as fol- 
lows from the BCS approximation. We note that at zero 
temperature the fraction of superconducting particles is 
high even at zero pairing strength this special point cor- 
responds to the absence of two-body interactions which 
results in pair- wise Fermi occupation of time-reversed or- 
bitals. 
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Figure 3: (Color online) The contour plot of the order param- 
eter as a function of the pairing strength and temperature. 
Half occupied 12-level system is shown. The line separates 
normal and paired regions based on the BCS equation. 

Throughout this work we mainly discuss systems with 
an even particle number; we found that the difference 
between odd and even systems in the critical region of 
interest is small. Most of the distinction occurs at zero 
temperature where degeneracy of an odd-particle ground 
state and non-zero spin are important. This can be seeing 
in Fig. |4]where we compare the entropy and specific heat 
as a function of temperature for N = 11 and N = 12 12- 
level ladder systems. 

The transition to the thermodynamic limit is explored 
for a two-level system in Fig. [H Unless noted otherwise, 
in our study we select exactly half-occupied systems with 
N = £1/2. The region of interest is identified by the peak 
in heat capacity seen in Fig. E](b). With the increased 
particle number this peak becomes sharper as expected in 
the macroscopic limit, where the phase transition is rep- 
resented by a discontinuity. Another interesting remark 
can be made about the location of the peak. Following 
Eq. (fTS")) within the BCS approximation the location of 
the phase transition for a half-occupied two-level model 
does not depend on size of the valence space O, at V = 1 
the BCS prediction is T c — 0.455. As seen from the figure 
this is not exactly correct, for a small 10-particle system 
the peak appears at about T c = 0.35, and only with the 
increase in the particle number the peak moves right to 
the BCS predicted value, thus confirming the BCS as an 
exact theory in the macroscopic limit. 

In recent years analysis of poles in the complex temper- 
ature plane and the evolution of branches of these poles 
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Figure 4: (Color online) (a) Entropy and (b) Specific heat 
as a function of temperature, for an odd and even number of 
particles and V = 1. 
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Figure 5: (Color online) (a) Energy and (b) Specific heat as 
a function of temperature, for V = 1 and various number of 
particles N=10, 30, 50, and 100. 
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Figure 6: (Color online) The lowest branch of zeros computed 
for N=100 and G=1.00 is schematically shown. 



parison with the thermal BCS for a two-level model can 
be found in Ref. 

In what follows we use the classification of phase tran- 
sitions developed by Bormann et.al. [ll|. We introduce 
complex temperature as B = (3 + ir and numerically 
seek a set of zeros B% in the canonical partition function 
Z(Bi, N) = 0, since the function is real the zeros appear 
in complex conjugate pairs and we can limit the region 
of consideration to t > 0. The product expansion of the 
partition function in terms of zeros using the Weierstrass 
theorem gives 



z(B) = nH(i 



B_ 

B~i 



1 - 



(29) 



The DOZ in the complex temperature plane for the two- 
level system is shown schematically in Fig El 

The sets of zeros form branches (l8l. l2a|. in FigElonly 
the branch lowest to the real axis is shown. The size 
of the system determines the distance between neighbor- 
ing zeros which in macroscopic limit becomes continuous. 
Phase transitions are associated with branches crossing 
the real axis. Indeed the zeros in the partition function 
appear as poles in thermodynamic variables; for energy 
or heat capacity we have from lj29j) 



(E(B))=J2 



1 



1 



B t -(3 Bf-P 



(30) 



has attracted a lot of attention as a study and classifi- 
cation tool for mesoscopic phase transitions. The theory 
related to the distribution of zeros (DOZ) in fugacity 
of t he g rand canonical ensemble dates back to Yang-Lee 
[13, Hf. Later works [H, US HI extended it the to the 
complex temperature plane of the canonical ensemble. 
The method of classifications of mesoscopic phase tran- 
sitions, recently suggested in Ref. is based on the 
distribution of zeros near the real axis. Some of the inter- 
esting questions such as whether the nature of the phase 
transition changes as a function of size have been stud- 
ied with this approach. The first steps in the analysis of 
mesoscopic systems undergoing pairing phase transitions 
were done in Ref. @, the evolution of DOZ and com- 



In general, although there are no poles at the real axis, 
the derivative d k (In Z) j 'd(3 k ~ Ylj(&i ~ P)~ k ma y result 
in a divergent sum. As suggested in (lll | the classifica- 
tion of phase transitions in the Ehrenfest sense can be 
extended to a smaller system by considering how the dis- 
crete roots of the phase transition branch approach the 
real axis. By labeling the roots in the phase transition 
branch starting from the closest one to real axis, see Fig. 
El the crossing angle can be given as 

v = arctan . 

T2 - Tl 



8 



The power law that expresses the congestion of roots as 
they approach real axis at r — > determines the second 
parameter a as |<Bj+i — 2?i| ~ r r"- 

The first order phase transition, which in thermody- 
namic limit appears as a discontinuity in the first deriva- 
tive of the free energy corresponds to a vertical uniform 
approach of poles v = 0, a = 0. In other cases the 
transition is of the second order for < a < 1 or of 
a higher order if a > 1. This classification establishes a 
condition at which poles in sums of the form IpO)) and 
(|3l"T) accumulate a logarithmically divergent series. For a 
vertical approach, v = at the critical temperature the 
Pcr\ ~ j 1 /( Q+1 ) therefore fe-th derivative of the parti- 
tion function would lead to a divergent series if k > a + 1. 

To find poles in the complex plane we developed a nu- 
merical technique that uses analiticity of the above func- 
tions. We first determine the number of roots in a desired 
region using a contour integral 



(32) 



The line integration is fast and is done avoiding paths 
that go directly over the roots, this assures numerical 
stability and the real and integer result of Eq. Il32[) guar- 
antees the accuracy. Once the number of roots is known 
we use a method in the spirit of the Laguerre's polyno- 
mial root finding technique [63] . The problem is math- 
ematically analogous to the two-dimensional problem of 
electrostatics. In the numerical method we converge to 
a given "charge" in the presence of the field from other 
"charges" which is modeled via multipole expansion using 
the analytically known derivatives of the "field strength". 
The found roots are sequentially removed, namely bal- 
anced by the "charge" of an opposite sign. Depending on 
the starting point and the density of roots, the numeri- 
cal cancellation is not always perfect, and the same root 
may appear several times. Given that the total number 
of roots is known this problem is easily fixed by choosing 
a different starting point or by exploring a smaller region. 
In the calculations we stabilize the sum in the partition 
function by selecting scaling so that the largest term in 
the sum f22"]) equals unity. 

A series of plots where evolution of poles in the com- 
plex temperature plane as a function of the pairing 
strength is shown in Fig. [Jj The behavior of the heat ca- 
pacity as a function of temperature for each case is shown 
below to highlight the phase transition point. With no 
pairing, V = 0, the zeros are distributed along the two 
(almost) horizontal lines. Similar picture is seen at the 
pairing strength significantly below critical (V cr = 0.5 
at zero temperature from BCS). At about the critical 
strength, V — 0.4, a noticeable bifurcation occurs with 
the lower branch evolving toward the real axis. As pair- 
ing strength increases, the branches move down and more 
branches becomes visible in our figures; in Fig. [jjwe use 
the same temperature scale for all values of V. The lowest 
branch that approaches the real axis is associated with 
the phase transition. The latter is confirmed by the peak 



in the heat capacity that becomes sharper in cases with 
stronger pairing. 

In Fig. [8] the dependence of the critical parameters v 
and a on the pairing strength is addressed. Below the 
critical pairing strength the curves fluctuate, here, there 
is no phase transition and v and a can not be interpreted 
as critical parameters. At a sufficiently strong pairing 
interaction, however, the behavior of the parameters sta- 
bilizes showing a second order phase transition. 



B. Grand Canonical ensemble 

The grand canonical ensemble is of a special impor- 
tance in statistical mechanics, since the partition func- 
tion for non-interacting particles, Zo(/3, /Lt), is given by 
an analytical expression. The grand canonical partition 
function can be determined using the canonical one, 



Z(J3, = z N Z(/3, N), 



(33) 



N=0 



where fugacity z = exp(/3^i) is introduced. For non- 
interacting Fermi particles 

Zo(/3, V) =H(l + zexp(-Pe 1 )), 



where e\ is a single-particle spectrum. The above expres- 
sion results, for example, in the commonly used form for 
the occupation numbers 

rii = (1 + exp[/3(ei - • 

The grand canonical approach and the use of the above 
Fermi distribution for small systems with a fixed num- 
ber of particles is common, however may present serious 
problems. On the other hand, even for non-interacting 
particles computation of the microcanonical or canoni- 
cal partition function is difficult Investigation of 
the mesoscopic limit where statistical ensembles may no 
longer be equivalent is of interest here. 

The grand canonical ensemble is advantageous even 
when it comes to interacting systems, the partition func- 
tion can be expressed via diagrammatic summation. In 
relation to pairing we mention here a method first pro- 
posed in Ref. [63|, a more detailed discussion can be 
found in [(3(|. The full partition function for constant (or 
factorizable) pairing interaction can be obtained as 



poo 

Z = Z Q dt exp(-y(i)), 
Jo 

where function Y(t) is 



(34) 



y(t) = <-4£> 



1>0 



cosh 
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Figure 7: (Color online) Evolution of DOZ and C v in the complex temperature plane B — (3 + ir for 7V=100 particles in the 
half-occupied two-level system. There are number of poles near and exactly on the imaginary axis, they are of no interest to 
our discussion and are not shown. The poles for other systems are discussed in [2|] and references therein, the interpretation 
and the nature of branches is discussed in [fil], further in-depth exploration of the above model can be found in Ref. [25l |. 




Pairing strength 

Figure 8: (Color online) Parameter of the phase transition for 
the two level system with iV=100 particles as a function of 
pairing strength (a) a vs V and (b) transition angle v vs V. 



The most straightforward saddle point approximation to 
the integral |34|) leads to a saddle point t s , which we ex- 
press in terms of a gap parameter as A 2 = Gt s /(2(3). 
Thus, the saddle-point equation becomes a familiar gap 



equation of thermal BCS (fl4l) . and the thermal BCS 
theory represents the lowest order approximation of the 
grand canonical expression in Eq. f34|) . 

Various thermodynamic properties of the ladder sys- 
tem with 12 levels and 12 particles obtained with the ex- 
act calculation of the grand canonical partition function 
are shown in Fig. [9l the figure also includes compar- 
isons with the corresponding curves from the canonical 
ensemble, where applicable. The fluctuations in the par- 
ticle number as a function of temperature are shown in 
Fig. [9](a) . The value of this quantity levels at about two 
particles, a similar uncertainty on a level of one pair is 
present in the BCS theory. The particle uncertainty rela- 
tive to the system size ~ 2/N can be used to estimate the 
quality of the grand canonical ensemble in applications 
to particle-conserving mesoscopic systems. The results 
of comparisons between canonical and grand canonical 
ensembles for the entropy, excitation energy, and spe- 
cific heat as a function of temperature are shown in Fig. 
[9jb)-(d). The difference is quite small, and is consistent 
with the level of error from the particle non-conservation. 
Further comparison is shown in Fig. [10] where the en- 
ergy difference between canonical and grand canonical 
ensembles is plotted as a function of temperature. In the 
picket-fence model the discrepancy is noticeable, however 
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Figure 9: (Color online) Thermodynamic properties of the 
ladder system with 12 levels, 12 particles, and 1^=1.00 are 
shown as a function of temperature. The grand canonical 
ensemble is compared to the canonical, (a) Fluctuation in 
the number of particles in the grand canonical ensemble; (b) 
entropy; (c) specific heat; (d) excitation energy. 
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Figure 10: (Color online) The excitation energy difference 
between canonical and grand canonical statistical ensembles 
is shown as a function of temperature. Left panel corresponds 
to a ladder system and two-level model is on the right. 



it becomes relatively small in the two-level case with a 
much larger number of particles. The difference peaks ex- 
actly at the temperature of the phase transition (in both 
cases V = 1) where fluctuations are large. As seen in 
the two-level model for a large number of particles this 
region becomes narrow. Although for a ladder system 
the difference grows in the absolute scale, this behavior 
is associated with the extreme smallness of the system 
and the discrepancy per particle is still going to zero. 

The zeros of the analytic continuation of the grand 
partition function into the complex plain of chemical po- 
tential are of a separate interest. We start by defining 
these points with the set of complex numbers fii that for 



a certain temperature satisfy the equation Z(/3, fii) = 0. 
There are some features to be stressed here. The num- 
ber of principal roots is equal to the capacity of the 
fermion space f2. The grand canonical partition func- 
tion (|33|) is an Sl-th order polynomial in fugacity which 
leads to roots in the chemical potential that can be 
found with the standard numerical techniques for poly- 
nomials. The methods discussed in the context of the 
canonical partition function are also useful in this case. 
As the size of the system grows the roots increase in 
number and may form branches that can approach the 
real axis. This describes the mesoscopic phase transi- 
tion within the Yang-Lee picture [27t l28j . The accumula- 
tion of roots near the real axis, similarly to the canonical 
ensemble discussed earlier, represents a phase transition 
marked by the discontinuity in a certain order derivative 
of the grand canonical partition function with respect to 
chemical potential. This leads to the discontinuity in the 
pressure- volume diagram [67] and in the thermodynamic 
potential as a function of the particle number, namely 
condensation. Based on the well known properties of 
the Bose gas the appearance of such third order transi- 
tion could be a good evidence for the Bose-Einstein 
pair condensation. Whether with the increased pairing 
strength or in a certain limit of temperature the Cooper 
pairs become dynamically equivalent to bosons and form 
a condensate and if there is a crossover region is an in- 
teresting and important question (H5 , l69| . 

Before addressing the results of this study we discuss 
some of the expected features that can be inferred from 
the partition function (|34| . Within the BCS approxima- 
tion the integral 11341) is given by the single saddle point 
value 



2/3A 2 , T-r 

Zscs = Z exp( — ) _[_[ 

1>0 



G 



cosh 



(««) 



cosh 



The Yang-Lee zeros of the above expression are zeros of 
the hyperbolic cosine and for each single particle energy 
ei an infinite series of roots labeled by integer n can be 
obtained 



2n+ 1 



(35) 



The evolution of DOZ of grand canonical ensemble in 
the complex plane of chemical potential for a fixed pairing 
strength G and various temperatures is shown in Figs. 
[TT1fT3l In all of the plots only the principal branch of 
roots with n = 0, Eq. (f35| . which is closest to the real 
axis is shown. 

In Fig. Ql] a somewhat high pairing strength V = 1 
is selected so that at low temperature the system is well 
in the superconducting phase. The resulting zeros are lo- 
cated along the horizontal line consistently with Eq. (|35| . 
As the temperature increases the two lines of roots move 
apart deeper into the complex plane; this trend is again 
in agreement with (|35|) . however the overall behavior of 
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the roots is no longer regular. The critical temperature in 
this system (from the peak in heat capacity), is T c = 0.38 
which coincides with a region where the behavior of DOZ 
changes. As seen from the figure there are no branches 
of substantial significance that cross real axis, indicating 
no phase transition. 

The following Figs. [32] and [13] repeat the same study 
with weaker and stronger pairing. The findings are simi- 
lar: at about critical temperature the DOZ changes from 
the two-line distribution, reflecting the BCS limit, to a 
more scattered set of roots moving away from the real 
axis at high temperature. 

This exact calculation is consistent with the similar 
study [69]. At temperatures below critical and with 
strong pairing, Fig. [13], there are small symmetric 
branches of zeros that are directed toward the real axis 
(although never reach it), they do not appear to result 
in any transitional behavior and their significance is un- 
clear. Our models lack the explicit spatial degree of free- 
dom and it is likely that the BCS-BEC transition that 
reflects the change in the physical size of the Cooper pairs 
is simply impossible here. 
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Figure 11: (Color online) The distribution of zeros (DOZ) of 
grand canonical ensemble in the complex chemical potential 
plane for a two levels system with V=50 and V=1.00 for var- 
ious temperatures indicated. The T c = 0.38 for this system. 



The grand canonical partition function f34|) is use- 
ful for understanding DOZ in the complex temperature 
plane, although canonical and grand canonical ensembles 
are not fully equivalent. Balian and Langer [7(j have de- 
termined that the zeros approach the real axis at an angle 
v = 7r/4 and their density tends linearly to zero, a = 1, 
showing a second order transition. 
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Figure 12: (Color online) Same as Fig. 1111 except V = 0.5 
and corresponding T c — 0.17. 
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Figure 13: (Color online) Same as Fig. [TT] except V = 2 and 
thus T c = 0.8. 



C. Microcanonical ensemble 

The microcanonical ensemble is often assumed to be 
the most physically appropriate statistical treatment of 
a closed system. There is a number of theoretical works 
as well as direct experimental studies of nuclear thermo- 
dynamics in the microcanonical ensemble [l9L l6ll. fTll. [7"3| . 
The density of states (DOS) p is the primary element in 
the approach. Regrettably, the formal definition given 
earlier fTS]! is not appropriate per se, the density of states 
as well as the normalization in the discrete spectrum re- 
quires some averaging energy interval. For most of this 
study we chose not to implement a traditional binning 
method substituting it by the propagator-type approach, 
where an artificially inserted small width a (the same for 
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all states) results in the Lorentzian-type smoothing of 
every peak. The derivatives of the DOS are then calcu- 
lated based on the analytic derivatives of the Lorentzian 
which provides an additional stability. With this proce- 
dure the DOS p{E) is obtained. The averaging width 
a is an artificial parameter that introduces thermody- 
namic averaging, the results may strongly depend on this 
parameter when it is smaller than the average level spac- 
ing. This parameter is not necessarily a disadvantage, 
to the contrary, it allows us to zoom at the energy scale 
of interest. Within this work we select a — 0.5 — 1.0 in 
single-particle level spacing units. This is most reason- 
able micro scale and can be compared with the resolution 
scale of the canonical ensemble where energy fluctuations 
are at about 10. The a interval versus level spacing can 
be interpreted as the number of states needed to obtain a 
statistically equilibrated value for observable quantities, 
in the limit of quantum chaos a single state is sufficient 
[iiL [73| , on the other hand as discussed below pure pair- 
ing due to seniority conservation is poorly equilibrated 
and many states must be included. The latter fact influ- 
enced our choice of a. 

The entropy in the microcanonical ensemble is 
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S(E) = hi p(E), 
and the temperature can be defined as 



(36) 



(37) 



which does not depend on the normalization that is used 
for the DOS. 

In Fig. [I4]the temperature is shown as a function of ex- 
citation energy for all three ensembles. The microcanon- 
ical curve with a = 1 shows several low-lying peaks that 
can be identified with the pair breaking 0, El III • The 
seniority is a conserved quantum number for pure pairing 
interaction, nevertheless these peaks survive in the pres- 
ence of all interactions as was shown in Ref. |4Sil7lll74l]. 
The corresponding oscillations in the heat capacity and 
especially the regions where this quantity is negative can 
be associated with the paired to normal phase transi- 
tion which takes place in the pair-by-pair microscopically 
fragmented process. The canonical and grand canonical 
ensembles due to thermalization created by the heat bath 
smooth out these microscopic features into a single phase 
transition, and importantly, the same can be done in the 
microcanonical treatment by choosing a large averaging 
interval a. Already at a — 5, which is about half the en- 
ergy fluctuation in the canonical ensemble, the peaks in 
Fig. Q3] disappear and microcanonical approach becomes 
similar to canonical and grand canonical. 

The macroscopic limit of the microcanonical ensemble 
is considered in Fig. [15] where entropy as a function of 
excitation energy in the two-level system for various N is 
shown. The comparison of microcanonical, canonical and 
grand canonical treatments indicates that they become 
identical in the thermodynamic limit. The discrepancy 
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Figure 14: (Color online) Temperature as a function of en- 
ergy in three different statistical treatments: canonical, grand 
canonical and microcanonical. The ladder system with 12 lev- 
els, 12 particles and V=1.00 is used for this study. The results 
for the microcanonical ensemble are plotted with two different 
choices of energy window a = 1 and 5. 
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Figure 15: (Color online) Comparison of entropy as a function 
of excitation energy for the two-level system in three ther- 
modynamic ensembles (a = 0.5 in microcanonical). Pairing 
strength V=1.00. (a) For N=20 particles; (b) N=50 particles; 
(c) N=100 particles. 



at high energy is related to the finite space where in the 
microcanonical case the density of states becomes zero 
once the energy exceeds the maximum possible value for 
the model space. The model space limitation is a natural 
cut-off for all ensembles at high energy. 

In contrast to the canonical and grand canonical en- 
sembles where thermalization is provided by the external 
heat bath the thermalization is a serious question in the 
microcanonical treatment ¥7a, I76JI. It has been noted in 
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Refs. [20, \5M, [ZJ| that pairing interactions do not pro- 
vide sufficient thermalization. Particle-particle interac- 
tions of the pairing type only are not sufficient to fully 
mix states and thermodynamically equilibrate the sys- 
tem. Temperature determined microscopically ((371 is 
inconsistent with the one that comes from the occupa- 
tion numbers of individual single-particle levels see Ref. 
[2(J . 1 5 J . l74l |. This property of pairing makes the micro- 
canonical treatment special. The question of thermaliza- 
tion in systems with pure pairing is rather academic; as 
it has been shown in [2d] and references therein, at an ar- 
bitrary weak non-pairing interactions the equilibration is 
immediately restored. The significant role of non-pairing 
interaction was explored in recent work The mag- 
netic field discussed below can also provide the needed 
thermalizing effect. The sharp differences in the sta- 
tistical approaches seen in Fig. [14] suggest to look for 
an alternative treatment and tracking of the transitional 
behavior which would not introduce the heat bath, en- 
ergy averaging, or the particle number uncertainty but at 
the same time is statistically equilibrated. The Invariant 
Correlational Entropy in the next section provides a tool 
that satisfies this criteria. 



V. INVARIANT CORRELATIONAL ENTROPY 

The Invariant Correlational entropy (ICE) [H3, is a 
powerful method that allows phase transition features to 
be explored on a quantum mechanical level. Expanding 
the formal definitions of Sec. IIV[ the ICE for an indi- 
vidual eigenstate a is computed by averaging the density 
matrix over the interval of pairing strength 



I a = -TrQ^lnp"), p% k , = (k\a) (a\k>) , 

here k is a basis state. The final result due to the trace 
operation is basis independent. In Fig. [16] we show the 
invariant correlational entropy for all states in the paired 
N = 10 two-level system as a function of the excita- 
tion energy of the corresponding state. The ICE fluctu- 
ates from state to state and the curve shown has been 
smoothed. The enhancement of the ICE in the region 
between and 6 energy units of excitation signals a tran- 
sitional behavior. Indeed, a lower figure that shows the 
specific heat as a function of energy for the same system 
in the canonical ensemble reveals a consistent trend. The 
advantage of the ICE is that unlike canonical (or grand 
canonical) ensemble it needs no heat bath and maintains 
an exact particle conservation; on the other hand, it is 
not prone to equilibration and thermalization issues since 
those are established by the fluctuations of the pairing 
strength. 



VI. PAIRING AND MAGNETIC FIELD 

The presence of magnetic field is known to influence 
the physics of the pairing state and the pairing phase 
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Figure 16: (Color online) Invariant correlational entropy 
(ICE), upper panel; and the specific heat in the canonical 
ensemble, lower panel, are shown as a function of excitation 
energy for the two-level system with N=10 particles and V=l. 



transition. In this section we extend our study by show- 
ing the changes to the results brought by the presence of 
the field. The Hamiltonian to be considered here is 



Hb = H — gJB, 



(38) 



where J is the angular momentum of the state and B is 
the magnetic field. Without loss of generality we choose 
units of the magnetic field so that the gyromagnetic ratio 
g=l. The introduction of the magnetic field does not 
require a new diagonalization of the Hamiltonian. All 
eigenstates shift in energy according to their magnetic 
quantum number M, with the quantization axis along 
the field B. The problem outlined with the Hamiltonian 
([38]) is identical to the cranking model of rotating nuclei 
where H u = H — J ■ w with uj representing rotational 
frequency. Thus, a reader interested in rotations should 
understand "magnetic field" as a "rotational frequency". 

In the case of a single particle on one level the magne- 
tization is a textbook example 



(2j 



1 troth ( ~(2j 
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The spin fluctuation x(P, B) = (M 2 ) - (M) s 
spin susceptibility, is given by 



(39) 
related to 



X 



= — csch 



(|)_ w + 1) w(HH). 



(40) 

Here x = gB/T. The generalization of these results to 
the cases with many levels is straightforward. 

In our study below we assume that the degeneracy in 
the single-particle spectrum is due to the corresponding 
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value of the angular momentum j = w — 1/2 for the two- 
level model and j = 1/2 for all levels in the ladder sys- 
tem. For each seniority we deduce the number of states 
with certain angular momentum which in turn allows us 
to calculate statistical partition functions. The analyti- 
cally computed summations over the magnetic quantum 
numbers such as in Eqs. (|40|l and ([39]) speed up the cal- 
culations. 

The destruction of the superconducting state occurs 
because of the two somewhat related phenomena. A mag- 
netic field causes the breaking of superconducting pairs 
due to the lowering in energy of the spin aligned states. 
The estimate for the critical field in this case can be 
obtained by comparing the energy of the paired ground 
state with the seniority s — 2 aligned spin state with spin 
J. The latter is by 2A higher in energy at zero field and 
pairing becomes unfavorable when magnetic field exceeds 
B cr = 2A/(gJ). In our models, equations such as fTS]) or 
(fl5]) can be used for an estimate. The second reason is 
the change in the energy of the normal state reflecting the 
Pauli spin paramagnetism. In our models (half-occupied 
for the two-level case) the field above B cr = Ae/(2gj) 
will promote the particle-hole excitations across the gap 
between the single particle levels. It has been suggested 
7811 and experimentally confirmed, for example in Ref. 
79j, that due to these phenomena a sufficiently strong 
magnetic field could change the transition type from sec- 
ond to first order. The situation can be quite complex 
in mesoscopic systems where even without the field the 
classification can be somewhat difficult. 

The features of a 20-particle half-occupied two-level 
system are shown in Fig. [T7] as a function of tempera- 
ture for a set of different values of magnetic fields. Except 
for B = 0, for all curves in this figure the field exceeds 
both of the above critical values (B cr ^0.1). The behav- 
ior of the heat capacity illustrates the disappearance of 
the phase transition for all field strength shown, except 
for B = where a sharp peak is present. The average 
magnetization Fig. flTTc) is exactly zero in the absence of 
field and for high fields starts almost at the saturation. 
With increased temperature magnetization drops down. 

The set of fields below critical is shown in the next fig- 
ure (fT8l) . The behavior of the magnetic spin fluctuations 
is regular at B — 0, at weak fields this quantity exhibits a 
sharp peak at low temperatures, for strong fields the reg- 
ular behavior is again restored. The same peak is present 
in the spin susceptibility curve which is only by a factor 
of 1 different from X- The critical behavior is associated 
with the corresponding behavior of the magnetization, 
see Figs. [19] and [2Ql The peak in the heat capacity is 
reduced and shifted to lower temperature, showing that 
at non-zero field a superconductor has lower critical tem- 
perature. 

These results are quite robust. In Fig. [19] we show 
the same study repeated for the ladder model. The mag- 
netization (upper left panel in Fig. [19]) is zero for no 
field, it shows a peak in the region of the field strengths 
that corresponds to competition between the paired state 
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Figure 17: (Color online) Thermodynamics of the two level 
system with iV=20 and 1^=1.00 in the magnetic field. The (a) 
spin fluctuations x\ (b) specific heat; (c) magnetization (ratio 
to the maximum possible value); and (d) energy are shown as 
functions of temperature for different field strengths. For this 
model T cr = 0.46 at B = 0. 
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Figure 18: (Color online) The same system as in Fig. 1171 but 
concentrating on the low magnetic fields, (a) spin fluctua- 
tions, (b) specific heat. 



and magnetic orientation (B = 5 and 7 curves), and fi- 
nally at strong fields and low temperature it saturates to 
(M) = 6, the maximum alignment state. The spin fluctu- 
ations, shown on the left lower panel, at zero temperature 
and no field, are consistent with the typical results [80| . 
The sharp peak at magnetic fields below critical again re- 
flects the transition associated with magnetic alignment. 
The emergence of the peak that in thermodynamic limit 
would become a discontinuity in the otherwise continuous 
curve at B = can be interpreted as the change in the 
type of the phase transition. The evolution of the transi- 
tion point as a function of temperature is seen in the heat 
capacity curve which again shows lowering of the tran- 
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Figure 19: (Color online) Thermodynamics within the canonical ensemble of a 12-level ladder system with N = 12 particles in 
the magnetic field. The pairing strength is V = 1. The panels include plots of magnetization, entropy, magnetic susceptibility 
and specific heat as a function of temperature. The curves correspond to five different values of the magnetic field B=0, 5, 7, 
10 and 20. 



sition temperature with the increased field. The finding 
can be summarized as the existence of two critical values 
for the magnetic field first one that corresponds to the 
change in the nature of the normal-to-superconducting 
transition; and at the second, higher value of the field, 
the paired state is no longer supported. The change in 
the nature of the phase transition is best seen in the spin- 
susceptibility curve which at low fields has no peak and 
acquires a peak consistent with the peak in heat capacity 
at higher values of the magnetic field. We associate this 
behavior with the analogous situation in the macroscopic 
superconductor where the second order phase transition 
becomes first order in the non-zero magnetic field. 

Recently, sharp differences between the systems with 
odd and even particle numbers have been discussed in 
the literature [H3, We address this in Fig. [201 where 
we show the same study as in Fig. [Tjj]but with N = 11 
particles. The primary difference between these results 
is that the ground state is degenerate and both magneti- 
zation and entropy are non-zero at low temperature and 
low field. Otherwise the results are almost identical. We 
conducted similar calculations for a two-level model with 



N = 19 particles but decided not to show the results 
because the difference between N = 20 and N = 19 is 
almost impossible to notice (except for the entropy and 
magnetization at zero temperature and zero field). 

The presence of the external magnetic field certainly 
has an effect on the distribution of zeros in the complex 
temperature plane. This evolution is explored in Fig. [2T1 
where the motion of roots is traced as the magnetic field 
is increased in small increments. The initial distribution 
of roots is connected with a line marked by B — 0. The 
second line connects the roots at B — 0.01. The gradual 
rotation of the branch responsible for the phase transition 
is seen, which eventually, at high fields, no longer orients 
the roots toward the real axis. Based on a similar picture 
but for a large system with N = 100 particles that was 
studied earlier in Fig. [22j we show the change in the 
critical parameters associated with this motion in the 
presence of the field. Interestingly, both a and angle v are 
approaching zero which marks the change in the phase 
transition type from second to first order. Unfortunately, 
the zero a is not reached since the field strength becomes 
larger then the critical (here B cr ~ 0.01) and the paired 
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Figure 20: (Color online) Thermodynamics within the canonical ensemble of a 12-level ladder system with 11 particles in the 
magnetic field. The pairing strength is V = 1. The panels include plots of magnetization, entropy, magnetic susceptibility and 
specific heat as a function of temperature. The curves correspond to five different values of the magnetic field B=0, 5, 7, 10 
and 20. 



phase disappears. The plot in Fig. lj22j) is ended at this 
point since it is becomes impossible to identify a branch 
of roots relevant to the phase transition. 



Finally we mention a thermalizing role of the magnetic 
field which breaks high degeneracies of states; and as well 
as changing the phase transition globally it reduces the 
number of individual pair breaking transitions seen in 
the microcanonical treatment. In Fig. [23] the entropy 
in the microcanonical ensemble is plotted as a function 
of excitation energy for different strengths of the exter- 
nal magnetic field. The number of peaks associated with 
the pair breaking is ten at zero field and this number 
becomes smaller as the magnetic field gets stronger, sim- 
ply because of the pair alignment. Thermal equilibration 
and the equivalence of ensembles are seen in the following 
Fig. El In contrast to the B — case in Fig. HH the dif- 
ference between canonical and microcanonical ensembles 
disappears at the magnetic field near critical. 



VII. SUMMARY CONCLUSIONS AND 
OUTLOOK 

Using an exact solution to the pairing problem in a 
picket-fence and two-level systems we addressed differ- 
ent views on the pairing phase transition, pair breaking, 
thermalization, behavior in the magnetic field or, equiv- 
alent^, rotation in the framework of the cranking model. 
We present a comprehensive study analyzing paired sys- 
tems with various tools and methods ranging from the 
BCS treatment to different thermodynamic approaches, 
including invariant correlational entropy and zeros of the 
partition functions in the complex plane of temperature 
and chemical potential. 

We found the microcanonical, canonical, and grand 
canonical thermodynamic approaches to be different 
when applied to small systems, although as expected they 
are consistent in the macroscopic limit. The grand canon- 
ical and canonical treatment are surprisingly close to each 
other even in the cases with a relatively small particle 
number. We find the microcanonical treatment to be the 
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Figure 21: (Color online) The distribution of zeros without 
and with the presence of external magnetic field for two level 
system N=60 and G=1.00 . 



Figure 23: (Color online) Entropy in the microcanonical pic- 
ture is shown as a function the excitation energy for several 
values of the magnetic field below critical. The two-level sys- 
tem with 20 particles and V = 1 pairing strength was used. 
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Figure 22: The evolution of critical parameters a and v as a 
function of applied magnetic field B in system with iV=100 
particles and y=1.00 




Energy 



most adequate approach for closed small systems, it al- 
lows for both global and relatively local, in terms of the 
energy scale, consideration of the statistical properties. 
The averaging energy window needed for determination 
of the density of states gives a broad freedom to the mi- 
crocanonical approach. The corresponding energy fluctu- 
ation in canonical and grand canonical treatments is too 
large in small systems and smooths out many significant 
statistical features such as individual pair breaking ob- 
served in experiment [7l|]. In thermally equilibrated sys- 
tems the energy window can be as small as few times the 
level spacing since in the full quantum chaos an individ- 
ual state is a representative of the surrounding statistical 



Figure 24: (Color online) The same system as in Fig. [23] 
the temperature as a function of energy is compared in mi- 
crocanonical and canonical ensembles at the magnetic field 
strength B = 0.05. The difference at high energy is due to a 
finite model space. 

properties [73]. The idealization of interaction limiting 
them to pairing only represents a dangerous problem: 
the pairing forces exclusively cannot establish full equili- 
bration, this necessitates a larger thermodynamic energy 
window. 

The Invariant Correlational Entropy that relies on fluc- 
tuations in the pairing strength as a source of equilibra- 
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tion appears to be a powerful statistical tool, capable of 
exploring most of the features inherent separately to mi- 
crocanonical, canonical and grand canonical ensembles. 
This tool is particularly important in identifying phase 
transition regions in mesoscopic systems. 

Turning to the properties of the phase transition we 
found as mentioned above the microcanonical treatment 
to be distinct from other ensembles. In thermodynamic 
limit, however, all treatments agree. In the further study 
of phase transitions we discussed the distribution of zeros 
of the canonical partition function in the complex tem- 
perature plane. We developed and implemented a nu- 
merical method for counting and finding all of the com- 
plex roots in a given region. In agreement with the ear- 
lier findings 0, HH, we observe branches of roots and 
study the properties of the one that approaches the real 
axis. The behavior of the roots is consistent with the 
second order phase transition as classified in [H, HH and 
confirms similar macroscopic results [701 ] . 

The recent interest to the crossover region between su- 
perconducting pairing and Bose-Einstein condensation of 
pairs prompted us to consider the potential condensation 
by looking for zeros of chemical potential in the grand 
canonical partition function. The presence of such ze- 
ros near real axis would hint on the condensation phe- 
nomenon. We did not find significant branches of roots 
evolving toward the real axis, and no critical behavior 
was observed in thermodynamics. It is likely that our 
model with no explicit spatial degrees of freedom is not 
appropriate for these questions. 

The last chapter of our work is devoted to interesting 



study of the mesoscopic phase transition in the presence 
of magnetic field. It is fully equivalent to rotations within 
cranking model. We found that there is a resemblance be- 
tween observed mesoscopic properties and those known in 
the macroscopic physics of superconductors. At low field 
the normal and superconducting phases are separated by 
the second order phase transition. In the next region 
of higher magnetic field the normal and superconducting 
phases are separated by the transition of a different na- 
ture associated with a simultaneous peak in spin suscep- 
tibility end enhanced spin fluctuations. Finally, at even 
higher fields a superconducting state is not supported at 
all. We conjecture that this behavior is a mesoscopic 
manifestation of the second to first order change in the 
transition type known in the thermodynamic limit. We 
also traced the evolution of zeros in the canonical parti- 
tion function as a function of magnetic field. We found 
that the classification of transition type as suggested in 
Ref. [11] is consistent with the above argument. 
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